#include "Equation.h"

template <typename Mesh_Type, typename Element_Type>
class RealPossion1 : public FEM2D_Possion<Mesh_Type, Element_Type,2>,
                     public FEM2D_Dirichlet_Condition<Mesh_Type, Element_Type,2>,
                     public FEM2D_CG_Solver<Mesh_Type, Element_Type>,
                     public FEM2D_IrregularMesh<Solution, Mesh_Type, Element_Type, SpMat>
{
};

template <typename Mesh_Type, typename Element_Type>
class RealPossion2 : public FEM2D_Possion<Mesh_Type, Element_Type,2>,
                     public FEM2D_Dirichlet_Condition<Mesh_Type, Element_Type,2>,
                     public FEM2D_CG_Solver<Mesh_Type, Element_Type>,
                     public FEM2D_RegularMesh<Solution, Mesh_Type, Element_Type, SpMat>
{
};

double u(const Point<2, double> &_p)
{
    return std::sin(M_PI * _p[0]) * std::sin(2.0 * M_PI * _p[1]);
};

double f(const Point<2, double> &_p)
{
    return 5.0 * M_PI * M_PI * u(_p);
};

int main(int argc, char *argv[])
{

    RealPossion1<Easymesh, P1Element<double>> pro1;
    RealPossion1<Easymesh, P2Element<double>> pro2;
    RealPossion1<RegularMesh, P1Element<double>> pro3;  //用无格式网格输出方法输出格式化网格
    RealPossion1<RegularMeshP2, P2Element<double>> pro4;
    RealPossion2<RegularMesh, P1Element<double>> pro5;
    RealPossion2<RegularMeshP2, P2Element<double>> pro6;

    pro1.Initial("data/D");
    pro1.Take_one_step(1, f, &u);
    pro1.L2_error(4, u);
    pro1.Output("pro1");

    pro2.Initial("data/D");
    pro2.get_mesh().build_dofs();
    pro2.Take_steps(2, f, &u);
    pro2.L2_error(6, u);
    pro2.Output("pro2",4);

    int nx = 4, ny = 4;
    if (argc > 2)
	{
		nx = atoi(argv[1]);
		ny = atoi(argv[2]);
	}
    Point<2, double> x0y0({0, 0}), x1y1({1, 1});

    pro3.Initial();
    pro3.get_mesh().set_lbc(x0y0);
    pro3.get_mesh().set_ruc(x1y1);
    pro3.get_mesh().set_nxny(nx, ny);
    pro3.Take_one_step(1, f, &u);
    pro3.L2_error(4, u);
    pro3.Output("pro3");

    pro4.Initial();
    RegularMeshP2 &mesh = pro4.get_mesh();
    mesh.set_lbc(x0y0);
    mesh.set_ruc(x1y1);
    mesh.set_nxny(nx, ny);
    pro4.Take_one_step(2, f, &u);
    pro4.L2_error(6, u);
    pro4.Output("pro4",4);

    pro5.Initial(x0y0, x1y1, nx, ny);
    pro5.Take_one_step(1, f, &u);
    pro5.L2_error(4, u);
    pro5.Output("pro5");

    pro6.Initial();
    pro6.Initial_Regularmesh(x0y0, x1y1, nx, ny);
    pro6.Take_one_step(2, f, &u);
    pro6.L2_error(6, u);
    pro6.Output("pro6",2);

    return 0;
};